必要パッケージの読み込み

library(tidyverse)
library(marelac)
library(lubridate)
library(readxl)
library(showtext)
library(flextable)
library(marelac)
library(grid)
library(formattable)
library(DT) #DataTableパッケージ

小型回流水槽実験の目的

 キャノピー構造を保持した状態で、海藻の光合成速度は、受ける流速によってどのように応答するか検証する。これまでの研究では、海藻の一部分を用いた研究が主流であったが、キャノピー構造は海水の流動と相互作用し、海藻への栄養塩の供給量に影響を及ぼしている可能性がある。そのため、キャノピー構造を有する場合、海藻の一部分を用いた結果と異なる可能性がある。

 本研究では、
 (1)キャノピー構造がキャノピー内の流速にどのように影響しているのか?
 (2)キャノピー構造を有する場合、海藻の光合成速度は、流速に対してどのように応答するのか?
 (3)その応答の要因は、キャノピーによって影響を受けた海水の流動によるものなのか?

 を明らかにすることを目的とする。

実験方法

光合成実験

  • 小型回流水槽に暗幕をし、明条件1h, 暗条件2hとした。
  • 密度条件ごとに個体を変えた(流速条件では同じ)
  • 実験開始前に、海水のDOを窒素ガスでばっ気することで、90%近くまで下げた。
  • 光合成実験の際には、海水の空気中の酸素ガスの交換が発生する。そのため、ガス交換速度を流速条件ごとに算出し、補正した。

流速測定実験

  • 流速はドップラー流速計(サンプリング間隔:200Hz、測定時間:3分)を用いて測定した。
  • 測定位置は、藻場の直前、leading edgeから15 cm, 45 cm, 藻場直後とし、測定深度は、底面から6 cm, 藻場上部の境界(流速、密度によって異なる)、16 cmとした。

キャノピーの形状変化の測定

  • キャノピーの流速による形状変化は、キャノピーの横方向からの投影面積で表現した。
  • カメラで流速ごとに撮影し、imagejを用いて面積を算出した。

光量の測定

  • 光量は約800\(\mu mol \: photons \: m^{-2}\:s^{-1}\)の一定。
  • 藻場底面の流速、密度ごとの光量は、ペンダントロガーにて測定(lux)

溶存酸素濃度の補正と総光合成速度の算出

溶存酸素濃度の生データの確認

それぞれの条件ごとに、測定初期の乱れを除去するため900秒を捨てる。変動速度測定に使用するのは、開始位置から40分間とする。 TRUEが明条件、FALSEが暗条件を示す。

## # A tibble: 1,968 x 11
## # Groups:   Hz, Light, Density [48]
##       ID    Hz Light Elapsed_time Density    DO  Temp Barometer   Sal
##    <dbl> <dbl> <lgl>        <dbl> <fct>   <dbl> <dbl>     <dbl> <dbl>
##  1     3   0.5 FALSE            0 High     6.82  21.3      101.  34.6
##  2     3   0.5 FALSE            1 High     6.81  21.3      101.  34.6
##  3     3   0.5 FALSE            2 High     6.82  21.3      101.  34.6
##  4     3   0.5 FALSE            3 High     6.81  21.3      101.  34.6
##  5     3   0.5 FALSE            4 High     6.82  21.3      101.  34.6
##  6     3   0.5 FALSE            5 High     6.8   21.4      101.  34.6
##  7     3   0.5 FALSE            6 High     6.81  21.3      101.  34.7
##  8     3   0.5 FALSE            7 High     6.81  21.3      101.  34.6
##  9     3   0.5 FALSE            8 High     6.81  21.4      101.  34.6
## 10     3   0.5 FALSE            9 High     6.81  21.4      101.  34.6
## # ... with 1,958 more rows, and 2 more variables: WW_sum <dbl>,
## #   Seaweed_volume_m3 <dbl>


生データから溶存酸素の変化速度を算出

#溶存酸素濃度の変化速度を算出するための関数作成
Glmmodel = function(y, x){
  summary(glm(y ~ x, family = Gamma(link = "log")))
}                        

Slope = function(x){               
  coef(x)[2]
} 

Slope_error = function(x){               
  coef(x)[4]
} 

#溶存酸素変化速度を算出算出
DO_analysis = Trident %>%
  group_by(Hz, Light, WW_sum, Seaweed_volume_m3, Density) %>%
  nest() %>% 
  mutate(Model = map(data, ~ glm(DO ~ Elapsed_time, data = ., family = Gamma(link = "log")))) %>% 
  mutate(Summary = map(Model, ~summary(.))) %>% 
  mutate(Slope = map_dbl(Model, ~Slope(.)),
         Stderr = map_dbl(Model, ~Slope_error(summary(.))))

#湿重量あたりの溶存酸素変化速度を算出
DO_analysis_WW = DO_analysis %>% 
 mutate(Slope = Slope / WW_sum, Stderr = Stderr / WW_sum) 

湿重量あたりの溶存酸素濃度の変化速度の確認

TRUEが明条件、FALSEが暗条件を示す。
暗条件のデータは、実質呼吸速度を示す。図を見ると、流速が増加するにしたがって呼吸速度が正になっている。測定開始時は、窒素ばっ気することで溶存酸素濃度を90%程に落している。そのため、空気中から、溶存酸素が海水に溶け込んでいることになる。流速が増加するにしたがって、呼吸速度よりも、空気中からの溶け込み速度(ガス交換速度)が高くなったことが考えられる。

## # A tibble: 48 x 10
##       Hz Light WW_sum Seaweed_volume_… Density data  Model Summary    Slope
##    <dbl> <lgl>  <dbl>            <dbl> <fct>   <lis> <lis> <list>     <dbl>
##  1   0.5 FALSE   453.          0.00063 High    <tib… <S3:… <S3: s… -4.15e-7
##  2   0.5 TRUE    453.          0.00063 High    <tib… <S3:… <S3: s…  2.32e-6
##  3   1   FALSE   453.          0.00063 High    <tib… <S3:… <S3: s… -5.82e-7
##  4   1   TRUE    453.          0.00063 High    <tib… <S3:… <S3: s…  2.55e-6
##  5   2   FALSE   453.          0.00063 High    <tib… <S3:… <S3: s… -4.68e-7
##  6   2   TRUE    453.          0.00063 High    <tib… <S3:… <S3: s…  2.96e-6
##  7   4   FALSE   453.          0.00063 High    <tib… <S3:… <S3: s… -5.04e-7
##  8   4   TRUE    453.          0.00063 High    <tib… <S3:… <S3: s…  3.26e-6
##  9   6   FALSE   453.          0.00063 High    <tib… <S3:… <S3: s… -2.50e-7
## 10   6   TRUE    453.          0.00063 High    <tib… <S3:… <S3: s…  3.12e-6
## # ... with 38 more rows, and 1 more variable: Stderr <dbl>

溶存酸素変化速度の補正:海水と空気中のガス交換速度にて

 小型回流水槽で実験する際、測器の取り付けや実験の効率を考えて、観測部分水路の上部は開放した状態で実験を行った。そのため、流速が速くなるにつれて、海水-空気間の酸素ガスの交換が活発となり、溶存酸素濃度の変動に影響を与えると考えられる。そこで、海藻を入れていない状態で、窒素ばっ気により、溶存酸素濃度を低下させた海水の溶存酸素濃度の経時変化を測定し、流速条件ごとのガス交換速度を算出した。
 実験を行った流速条件は6, 0.5, 4, 8, 20 Hzとした。

補正に使用するデータの確認。

#大気-海水間の酸素フラックスの算出-------------------------------------------------------------
#光合成実験3で測定したデータを使用する。回流水槽内に海藻は入っていない。
Trident_flux = read_csv("../回流水槽光合成実験3/Modified_data/DO_air_water_Flux_data.csv") 
print(Trident_flux, n = 10)
ggplot(Trident_flux)+
  geom_point(aes(x = Elapsed_time, y = DO, color = "DO"))+
  geom_point(aes(x = Elapsed_time, y = DO_saturation, color = "DO_saturation"))+
  facet_wrap("Hz")

## # A tibble: 3,177 x 11
##    Time                ID    Barometer    DO  Temp   Sal
##    <dttm>              <chr>     <dbl> <dbl> <dbl> <dbl>
##  1 2018-05-06 06:15:27 calb…      101.  6.88  18.9  34.7
##  2 2018-05-06 06:16:27 calb…      101.  6.88  19    34.7
##  3 2018-05-06 06:17:27 calb…      101.  6.88  19    34.7
##  4 2018-05-06 06:18:27 calb…      101.  6.89  19    34.7
##  5 2018-05-06 06:19:27 calb…      101.  6.89  19    34.7
##  6 2018-05-06 06:20:27 calb…      101.  6.89  19    34.7
##  7 2018-05-06 06:21:27 calb…      101.  6.9   19    34.7
##  8 2018-05-06 06:22:27 calb…      101.  6.9   19    34.7
##  9 2018-05-06 06:23:27 calb…      101.  6.9   19    34.7
## 10 2018-05-06 06:24:27 calb…      101.  6.91  19    34.7
## # ... with 3,167 more rows, and 5 more variables: Date <dttm>, Hz <dbl>,
## #   DO_saturation <dbl>, DO_diff <dbl>, Elapsed_time <dbl>


 以上のデータを以下の式に当てはめる。  

 \(DO_{t}\:=\:C\:*\:exp(-k\:*\:t)\:+\:DO_{sat}\)
 \(DO_{t}\):経過時間tにおける溶存酸素濃度
 \(C\):t = 0における溶存酸素濃度と飽和溶存酸素濃度の差
 \(k\):溶存酸素濃度のガス交換による変化速度
 \(DO_{sat}\):経過時間tにおける飽和溶存酸素濃度
 

#データを用いてKの算出(kが負になると酸素の放出、正は吸収)-----------------------------
library(nlstools)
model = function (k , C, t, DO_saturation){
  C * exp(-k * t) + DO_saturation
}

k = Trident_flux %>%
  filter(!ID == "ECSER09") %>% 
  group_by(Hz) %>% 
  nest() %>% 
  mutate(nlsout = map(data, function(x){nls(DO~model(k0,C0,Elapsed_time, DO_saturation), data = x, start = list(k0=0.11, C0=-4))}),
         k = map_dbl(nlsout, function(x){coef(x)[[1]]}))

#当てはまり確認
k
k$nlsout
## # A tibble: 5 x 4
##      Hz data                nlsout           k
##   <dbl> <list>              <list>       <dbl>
## 1  20   <tibble [241 × 10]> <S3: nls> 0.00554 
## 2   8   <tibble [395 × 10]> <S3: nls> 0.00213 
## 3   4   <tibble [665 × 10]> <S3: nls> 0.000864
## 4   0.5 <tibble [318 × 10]> <S3: nls> 0.000138
## 5   6   <tibble [620 × 10]> <S3: nls> 0.00244 
## [[1]]
## Nonlinear regression model
##   model: DO ~ model(k0, C0, Elapsed_time, DO_saturation)
##    data: x
##        k0        C0 
##  0.005539 -0.655192 
##  residual sum-of-squares: 0.006818
## 
## Number of iterations to convergence: 8 
## Achieved convergence tolerance: 0.00000138
## 
## [[2]]
## Nonlinear regression model
##   model: DO ~ model(k0, C0, Elapsed_time, DO_saturation)
##    data: x
##        k0        C0 
##  0.002126 -0.732538 
##  residual sum-of-squares: 0.01068
## 
## Number of iterations to convergence: 7 
## Achieved convergence tolerance: 0.000008053
## 
## [[3]]
## Nonlinear regression model
##   model: DO ~ model(k0, C0, Elapsed_time, DO_saturation)
##    data: x
##         k0         C0 
##  0.0008645 -0.7389470 
##  residual sum-of-squares: 0.01915
## 
## Number of iterations to convergence: 8 
## Achieved convergence tolerance: 0.00000008483
## 
## [[4]]
## Nonlinear regression model
##   model: DO ~ model(k0, C0, Elapsed_time, DO_saturation)
##    data: x
##         k0         C0 
##  0.0001376 -2.4797074 
##  residual sum-of-squares: 0.009746
## 
## Number of iterations to convergence: 7 
## Achieved convergence tolerance: 0.000009302
## 
## [[5]]
## Nonlinear regression model
##   model: DO ~ model(k0, C0, Elapsed_time, DO_saturation)
##    data: x
##        k0        C0 
##  0.002439 -4.372584 
##  residual sum-of-squares: 1.053
## 
## Number of iterations to convergence: 6 
## Achieved convergence tolerance: 0.00000157


すべてのHzのおいてkを算出していないため、Hzとkは直線回帰で近似できると仮定して、kを補足する

#算出した酸素フラックス速度定数kを持ちいて、実験中の酸素フラックスを見積もる:使用するモデルはパターン2のDO_saturation を変動モデル-----------------------------------------------------------------------------
#現在(2018-03-21現在)ではすべてのHzのおいてkを算出していないため、Hzとkは直線回帰で近似できると仮定して、kを補足する

Hz = data.frame(Hz = c(0.5, 1, 2, 4, 6, 8, 10, 15, 20))
k = k %>%
  full_join(Hz) %>% 
  arrange(Hz)

#Hz6の推定が悪いので、NAにする
k = k %>% mutate(k = ifelse(Hz == 6, NA, k))

#一次式で近似
k_pre = data_frame(Hz = k$Hz, k = predict(lm(k ~ Hz, k), data = data_frame(Hz = k$Hz),newdata=data.frame(Hz = k$Hz)))

ggplot()+
  geom_point(aes(x = Hz, y = k), data = k)+
  geom_line(aes(x = Hz, y = k), data = k_pre)


 kを用いて溶存酸素濃度を補正する

#kを用いてDO濃度を補正
DO_offset_Flux = Trident %>% 
  mutate(DO_saturation = gas_O2sat(S = Sal, t = Temp),
         DO_diff = DO_saturation - DO) %>% 
  left_join(k_pre, by = "Hz") %>%
  mutate(DO_Flux_min = DO_diff * k) %>% 
  group_by(ID, Hz, Light, Density) %>% 
  mutate(DO_Flux_min_cum = cumsum(DO_Flux_min),
         DO_offset_Flux = DO - DO_Flux_min_cum) %>%
  return()

print(DO_offset_Flux, n = 10)

ggplot()+
  # geom_point(aes(x = Elapsed_time, y = DO_Flux_min, color = "Flux"), data = DO_offset_Flux)+
  geom_point(aes(x = Elapsed_time, y = DO_offset_Flux, color = "補正後"), data = DO_offset_Flux)+
  # geom_smooth(aes(x = Elapsed_time, y = DO_offset_Flux, color = "offset_line"), method = "glm", data = DO_offset_Flux)+
  # geom_point(aes(x = Elapsed_time, y = DO_saturation, color = "Saturation"), data = DO_offset_Flux)+
  geom_point(aes(x = Elapsed_time, y = DO, color = "補正前"), data = DO_offset_Flux)+
  # geom_smooth(aes(x = Elapsed_time, y = DO, color = "raw_line"), method = "glm", data = DO_offset_Flux)+
  # geom_point(aes(x = Elapsed_time, y = DO_diff, color = "diff"), data = DO_offset_Flux)+
  # geom_point(aes(x = Elapsed_time, y = DO_Flux_min_cum, color = "cum"), data = DO_offset_Flux)+
  # geom_line(aes(x = Elapsed_time, y = k, color = "k"), data = DO_offset_Flux)+
  facet_grid(Density ~ Light + Hz)

## # A tibble: 1,968 x 17
## # Groups:   ID, Hz, Light, Density [48]
##       ID    Hz Light Elapsed_time Density    DO  Temp Barometer   Sal
##    <dbl> <dbl> <lgl>        <dbl> <fct>   <dbl> <dbl>     <dbl> <dbl>
##  1     3   0.5 FALSE            0 High     6.82  21.3      101.  34.6
##  2     3   0.5 FALSE            1 High     6.81  21.3      101.  34.6
##  3     3   0.5 FALSE            2 High     6.82  21.3      101.  34.6
##  4     3   0.5 FALSE            3 High     6.81  21.3      101.  34.6
##  5     3   0.5 FALSE            4 High     6.82  21.3      101.  34.6
##  6     3   0.5 FALSE            5 High     6.8   21.4      101.  34.6
##  7     3   0.5 FALSE            6 High     6.81  21.3      101.  34.7
##  8     3   0.5 FALSE            7 High     6.81  21.3      101.  34.6
##  9     3   0.5 FALSE            8 High     6.81  21.4      101.  34.6
## 10     3   0.5 FALSE            9 High     6.81  21.4      101.  34.6
## # ... with 1,958 more rows, and 8 more variables: WW_sum <dbl>,
## #   Seaweed_volume_m3 <dbl>, DO_saturation <dbl>, DO_diff <dbl>, k <dbl>,
## #   DO_Flux_min <dbl>, DO_Flux_min_cum <dbl>, DO_offset_Flux <dbl>


補正した溶存酸素濃度を用いて、変化速度を再計算

DO_analysis_k = DO_offset_Flux %>%
  group_by(Hz, Light, WW_sum, Seaweed_volume_m3, Density) %>%
  nest() %>% 
  mutate(Model = map(data, ~ glm(DO_offset_Flux~Elapsed_time, data = ., family = Gamma(link = "log")))) %>% 
  mutate(Summary = map(Model, ~summary(.))) %>% 
  mutate(Slope = map_dbl(Model, ~Slope(.)),
         Stderr = map_dbl(Model, ~Slope_error(summary(.))),
         Temp_mean = map_dbl(data, function(x){
           x %>%
             summarise(Temp = mean(Temp)) %>% 
             pull(Temp)}),
         Temp_sd = map_dbl(data, function(x){
           x %>%
             summarise(Temp = sd(Temp)) %>% 
             pull(Temp)}))

#湿重量あたりの溶存酸素変化速度を算出
DO_analysis_k_WW = DO_analysis_k %>% 
 mutate(Slope = Slope / WW_sum, Stderr = Stderr / WW_sum) 


補正後の溶存酸素変化量の確認。
呼吸速度は、0以下の値を取るように補正された。光合成速度は山形の変動を示した。

## # A tibble: 48 x 12
##       Hz Light WW_sum Seaweed_volume_… Density data  Model Summary    Slope
##    <dbl> <lgl>  <dbl>            <dbl> <fct>   <lis> <lis> <list>     <dbl>
##  1   0.5 FALSE   453.          0.00063 High    <tib… <S3:… <S3: s… -4.17e-7
##  2   0.5 TRUE    453.          0.00063 High    <tib… <S3:… <S3: s…  2.32e-6
##  3   1   FALSE   453.          0.00063 High    <tib… <S3:… <S3: s… -6.09e-7
##  4   1   TRUE    453.          0.00063 High    <tib… <S3:… <S3: s…  2.52e-6
##  5   2   FALSE   453.          0.00063 High    <tib… <S3:… <S3: s… -5.47e-7
##  6   2   TRUE    453.          0.00063 High    <tib… <S3:… <S3: s…  2.88e-6
##  7   4   FALSE   453.          0.00063 High    <tib… <S3:… <S3: s… -6.37e-7
##  8   4   TRUE    453.          0.00063 High    <tib… <S3:… <S3: s…  3.20e-6
##  9   6   FALSE   453.          0.00063 High    <tib… <S3:… <S3: s… -4.34e-7
## 10   6   TRUE    453.          0.00063 High    <tib… <S3:… <S3: s…  3.07e-6
## # ... with 38 more rows, and 3 more variables: Stderr <dbl>,
## #   Temp_mean <dbl>, Temp_sd <dbl>

溶存酸素変化速度の補正:水温による補正

本実験では、水温の制御なしで実験を行った。光条件、密度条件ごとの平均水温、及び分散は以下のようになった。
Light Density Temp_mean Temp_sd
FALSE Low 20.926 0.219
FALSE Middle 21.084 0.249
FALSE High 21.612 0.283
TRUE Low 21.408 0.214
TRUE Middle 21.465 0.251
TRUE High 22.029 0.261


補正は以下のアレニウスの式を用いて行う。Durocher and Allen (2012)を参照。

\(r^{'}\:=\:r\:*\:exp(-Ea\:/\:k\:(\:1/K\:-\:1/K0\:)\)

 \(r^{'}\):補正後(水温がk0のとき)の値
 \(r\):補正前(それぞれの値の水温)の値
 \(k\):ボルツマン定数
 \(Ea\):活性化エネルギー(海藻の状態や実験条件によって異なると考えられる)

海藻の生理的な応答がアレニウスの式の関係に従うと仮定し、K0のときの水温に合わせて補正する。 まず、補正前のデータもちいてEaを推定する。補正する水温は20 °C とした。

Eaの値が、以上に大きいため、この値で補正するには不適切であると考えられる。

bolzmann = 8.617e-5 #ボルツマン定数
K0 = 20 +  273.15 #補正水温K0は全ての条件で統一する。

#光条件ごとにEaを推定し、推定したEaでSlopeを補正
DO_analysis_k_WW_Ea = DO_analysis_k_WW %>% 
  group_by(Light) %>%
  mutate(K0 = K0,
         K = Temp_mean + 273.15,
         invK = (1 / K0 * bolzmann - 1 / K * bolzmann), 
         Slope = ifelse(Light == T, Slope, Slope*(-1)),
         Slslope = scale(log(Slope), scale = F)) %>% 
  nest() %>% 
  mutate(Model = map(data, function(x){lm(Slslope ~ invK, data = x) %>% summary()}),
         Ea = map_dbl(Model, function(x){coef(x)[2]}),
         data_temp = map2(data, Ea, function(x, y){
           x %>%
             mutate(Slope = Slope * exp(-y * invK))})) 

DO_analysis_k_WW_Ea
DO_analysis_k_WW_Ea %>% 
  select(Light, data_temp) %>% 
  unnest() %>% 
  ggplot()+
  geom_point(aes(x = Hz, y = Slope, shape = "補正前", color = Temp_mean),
             size = 3,
             data = DO_analysis_k_WW %>% 
               mutate(Slope = ifelse(Light == T, Slope, Slope*(-1))))+
  geom_point(aes(x = Hz, y = Slope, shape = "補正後", color = Temp_mean),
             size = 3)+
  facet_grid(Light ~ Density)

#Light,Density 条件ごとにEaを推定した場合。
DO_analysis_k_WW_Ea = DO_analysis_k_WW %>% 
  group_by(Light, Density) %>%
  mutate(K0 = K0,
         K = Temp_mean + 273.15,
         invK = (1 / K0 * bolzmann - 1 / K * bolzmann), 
         Slope = ifelse(Light == T, Slope, Slope*(-1)),
         Slslope = scale(log(Slope), scale = F)) %>% 
  nest() %>% 
  mutate(Model = map(data, function(x){lm(Slslope ~ invK, data = x) %>% summary()}),
         Ea = map_dbl(Model, function(x){coef(x)[2]}),
         data_temp = map2(data, Ea, function(x, y){
           x %>%
             mutate(Slope = Slope * exp(-y * invK))})) 

DO_analysis_k_WW_Ea
DO_analysis_k_WW_Ea %>% 
  select(Light, Density, data_temp) %>% 
  unnest() %>% 
  ggplot()+
  geom_point(aes(x = Hz, y = Slope, shape = "補正前", color = Temp_mean),
             size = 3,
             data = DO_analysis_k_WW %>% 
               mutate(Slope = ifelse(Light == T, Slope, Slope*(-1))))+
  geom_point(aes(x = Hz, y = Slope, shape = "補正後", color = Temp_mean),
             size = 3)+
  facet_grid(Light ~ Density)

## # A tibble: 2 x 5
##   Light data               Model                     Ea data_temp         
##   <lgl> <list>             <list>                 <dbl> <list>            
## 1 FALSE <tibble [24 × 15]> <S3: summary.lm> -286625485. <tibble [24 × 15]>
## 2 TRUE  <tibble [24 × 15]> <S3: summary.lm>  219348551. <tibble [24 × 15]>
## # A tibble: 6 x 6
##   Light Density data            Model                   Ea data_temp      
##   <lgl> <fct>   <list>          <list>               <dbl> <list>         
## 1 FALSE High    <tibble [8 × 1… <S3: summary.l…    -2.21e8 <tibble [8 × 1…
## 2 TRUE  High    <tibble [8 × 1… <S3: summary.l…     2.37e8 <tibble [8 × 1…
## 3 FALSE Low     <tibble [8 × 1… <S3: summary.l…    -1.85e7 <tibble [8 × 1…
## 4 TRUE  Low     <tibble [8 × 1… <S3: summary.l…    -3.46e8 <tibble [8 × 1…
## 5 FALSE Middle  <tibble [8 × 1… <S3: summary.l…    -1.74e8 <tibble [8 × 1…
## 6 TRUE  Middle  <tibble [8 × 1… <S3: summary.l…     2.23e8 <tibble [8 × 1…


次にEaを0.65で固定して補正する。
補正前後で大きな変化はない。水温の変動から考えても妥当か。

#全ての条件でEaを固定する。
FixedEa = 0.65 #1000にしてもほどんど変化なし。

DO_analysis_k_WW_FixedEa = DO_analysis_k_WW %>% 
  group_by(Light) %>%
  mutate(K0 = K0,
         K = Temp_mean + 273.15,
         invK = (1 / K0 * bolzmann - 1 / K * bolzmann), 
         Slope = ifelse(Light == T, Slope, Slope*(-1))) %>% 
  nest() %>% 
  mutate(Ea = FixedEa,
         data_temp = map2(data, Ea, function(x, y){
           x %>%
             mutate(Slope = Slope * exp(-y * invK))})) 

DO_analysis_k_WW_FixedEa
DO_analysis_k_WW_FixedEa %>% 
  select(Light, data_temp) %>% 
  unnest() %>% 
  ggplot()+
  geom_point(aes(x = Hz, y = Slope, shape = "補正前", color = Temp_mean),
             size = 3,
             data = DO_analysis_k_WW %>% 
               mutate(Slope = ifelse(Light == T, Slope, Slope*(-1))))+
  geom_point(aes(x = Hz, y = Slope, shape = "補正後", color = Temp_mean),
             size = 3)+
  facet_grid(Light ~ Density)

## # A tibble: 2 x 4
##   Light data                  Ea data_temp         
##   <lgl> <list>             <dbl> <list>            
## 1 FALSE <tibble [24 × 14]>  0.65 <tibble [24 × 14]>
## 2 TRUE  <tibble [24 × 14]>  0.65 <tibble [24 × 14]>

水温補正は、実験間の水温のばらつきが最大で約1°C 程であることから、それほど必要ないかもしれない。推定したEaの絶対値が大きくなったのは、Eaを算出する際に、密度条件や個体差によるSlopeの違いが、推定に反映されたからではないか。
よって水温補正は、固定したEaを用いた補正を使用する。


総光合成速度の算出

総光合成速度は、純光合成速度と総光合成速度の和で算出される。

NP = DO_analysis_k_WW_FixedEa %>% 
  select(Light, data_temp) %>% 
  unnest() %>% 
  filter(Light == TRUE) %>% 
  mutate(type = "NP") %>% 
  select(type, Hz, Density, WW_sum, Temp_mean, Seaweed_volume_m3, Slope)

RP = DO_analysis_k_WW_FixedEa %>% 
  select(Light, data_temp) %>% 
  unnest() %>% 
  filter(Light == FALSE) %>% 
  mutate(type = "RP") %>% 
  select(type, Hz, Density, WW_sum, Temp_mean, Seaweed_volume_m3, Slope)

GP = NP %>% 
  full_join(RP, by = c("Hz", "Density", "WW_sum", "Seaweed_volume_m3")) %>% 
  mutate(Slope = Slope.x + Slope.y, 
         type = "GP") %>% 
  select(type, Hz, Density, WW_sum, Seaweed_volume_m3, Slope)

#データの結合
V = 200 #L 回流水槽の容量
Slope_DO = bind_rows(NP, RP, GP) %>% 
  mutate(Slope = Slope*V, #単位をmg/l/min/gww -> mg/min/gww
         Slope = Slope*1000) #単位をmg/min/gww -> μg/min/gww

head(Slope_DO)
## # A tibble: 6 x 7
##   type     Hz Density WW_sum Temp_mean Seaweed_volume_m3 Slope
##   <chr> <dbl> <fct>    <dbl>     <dbl>             <dbl> <dbl>
## 1 NP      0.5 High      453.      21.7           0.00063 0.464
## 2 NP      1   High      453.      22.0           0.00063 0.505
## 3 NP      2   High      453.      21.6           0.00063 0.575
## 4 NP      4   High      453.      22.3           0.00063 0.640
## 5 NP      6   High      453.      21.9           0.00063 0.614
## 6 NP      8   High      453.      22.2           0.00063 0.588

流速データの補正

流速データのノイズを除去する。Jesson et al., 2013参照 ### ノイズ除去のための関数作成

# Read packages ----------------------
library(lubridate)
library(zoo)

# Declare the function --------------------------
#xには流速の残差が代入される。差分法による微分値を算出する。
gradient = function(x){
  n = length(x)
  j = 2:(n-1)
  c(x[2] - x[1],
    sapply(j, function(j, z){(z[j+1] - z[j-1])*0.5}, z = x),
    x[n] - x[n-1])
}

#advdataには、流速の生データが代入される。残りの引数はそれぞれの閾値を示す。
#閾値は、論文より引用して決めるが、便宜上以下の値を使用。
#返り値は、閾値から外れたデータをNAに変化したもの。
snr_corr_filter = function(advdata, snr.threshold=2, corr.threshold=60) {  
  # A filter for the SNR and CORR.
  require(tidyverse)
  snr = advdata %>% select(dplyr::starts_with("snr")) %>% mutate_all(funs(na=. < snr.threshold))
  corr = advdata %>% select(dplyr::starts_with("corr")) %>% mutate_all(funs(na=. < corr.threshold))
  tmp = as.matrix(advdata %>% select(dplyr::contains("vel")))
  tmp[as.matrix(snr %>% select(dplyr::ends_with("na")))] =NA
  tmp[as.matrix(corr %>% select(dplyr::ends_with("na")))] =NA
  tmp
}

interpolate = function(value) {
  # Conducts a cubic spline interpolation of removed spikes.
  x = which(!is.na(value))
  y = value[x]
  xi = seq_len(length(value))
  spline(x,y, xout=xi, method="fmm")$y
}
ellipse_radius=function(ab, theta, alpha=0, x0=0, y0=0) {
  # Use the polar equation of the ellipse to determine the radius.
  x=x0+ab[1]*cos(theta)*cos(alpha)-ab[2]*sin(theta)*sin(alpha)
  y=y0+ab[1]*cos(theta)*sin(alpha)+ab[2]*sin(theta)*cos(alpha)
  radius=ab[1]*ab[2]/(sqrt(ab[1]^2*sin(theta-alpha)^2+ab[2]^2*cos(theta-alpha)^2))
  data.frame(x,y,radius)
}
distance_radius=function(x,y, ab, alpha=0) {
  # Determine the distance to the point and the distance (radius)
  # to the edge of the ellipse
  theta = atan(y/x)
  distance=sqrt(x^2+y^2)
  theta[x<0]=theta[x<0]+pi
  radius = ellipse_radius(ab, theta, alpha=alpha)
  data.frame(theta, distance, radius=radius$radius)
}

make_ellipse = function(ab,alpha=0,x0=0,y0=0,size=500) {
  theta=seq(0,2*pi,length=size)
  x=x0+ab[1]*cos(theta)*cos(alpha)-ab[2]*sin(theta)*sin(alpha)
  y=y0+ab[1]*cos(theta)*sin(alpha)+ab[2]*sin(theta)*cos(alpha)
  data.frame(x,y)
}

#残差と残差の微分値(一階、二階)の相関図は、楕円形となる。
#楕円作図用関数
plotdata=function(f,ft,ftt,f_ft, ft_ftt, f_ftt,f_ft_bad, ft_ftt_bad, f_ftt_bad, theta) {
  require(ggplot2)
  ell_f_ft = make_ellipse(f_ft,0)
  ell_ft_ftt = make_ellipse(ft_ftt,0)
  ell_f_ftt = make_ellipse(f_ftt,theta)
  plot1=ggplot(data.frame(f,ft)) +
    geom_point(aes(f,ft, color=f_ft_bad)) +
    geom_path(aes(x,y), data=ell_f_ft)
  plot2=ggplot(data.frame(ft,ftt)) +
    geom_point(aes(ft,ftt, color=ft_ftt_bad)) +
    geom_path(aes(x,y), data=ell_ft_ftt)
  plot3=ggplot(data.frame(f,ftt)) +
    geom_point(aes(f,ftt, color=f_ftt_bad)) +
    geom_path(aes(x,y), data=ell_f_ftt)
  gridExtra::grid.arrange(plot1,plot2,plot3, ncol=2)
}
MAD = function(x) {
  # See Wahl (2003) Discussion of "Despiking Acoustic Doppler Velocimeter Data"
  # by Derek G. Goring and Vladimir I. Nikora, Journal of Hydraulic Engineering, 129(6), 484-487.
  # and Goring and Nikora (2003) Closure to "Depiking Acoustic Doppler
  # Velocimeter Data" by Derek G. Goring and Vladimir I. Nikora,
  # Journal of Hydraulic Engineering, 129(6), 487-488.
  # This function calculates the median absolute deviation. The value 1.483
  # makes it similar to the standard deviation.
  sigma = median(abs(x - median(x, na.rm=T)), na.rm=T)
  k = 1.483 # for normal distribution
  k*sigma
}
despike=function(f, good_data, figures=TRUE) {
  ft = gradient(f)
  ftt = gradient(ft)
  theta = atan(sum(f*ftt)/sum(f^2))
  fsigma = MAD(f)
  ftsigma = MAD(ft)
  fttsigma = MAD(ftt)
  nobs = length(f)
  lambda = sqrt(2*log(nobs)) # Universal threshold
  # Set the minor and major axes for the ellipse
  f_ft = c(lambda*fsigma, lambda*ftsigma)
  f_ftt = c(lambda*fsigma, lambda*fttsigma)
  ft_ftt = c(lambda*ftsigma, lambda*fttsigma)
  # Determine the distance and radius of the ellipse for each point
  f_ft_tr = distance_radius(f, ft, f_ft)
  f_ftt_tr = distance_radius(f, ftt, f_ftt, alpha=theta)
  ft_ftt_tr = distance_radius(ft, ftt, ft_ftt)
  # Determine of the point is beyond the ellipse
  f_ft_bad = apply(f_ft_tr [, c("radius","distance")], 1, diff) > 0
  f_ftt_bad = apply(f_ftt_tr [, c("radius","distance")], 1, diff) > 0
  ft_ftt_bad = apply(ft_ftt_tr[, c("radius","distance")], 1, diff) > 0
  # Flag the spikes to remove
  makeNA = unique(sort(c(which(f_ft_bad), which(f_ftt_bad), which(ft_ftt_bad))))
  f_tmp=f
  makeNA = makeNA[!(makeNA %in% good_data)]
  f[(makeNA)] = NA
  f[1] = ifelse(is.na(f[1]), 0, f[1]) #初期値がNAの場合エラーを返すのでNAの場合は便宜上0を代入
  # f_new = interpolate(f)
  f_new = zoo::na.locf(f,na.rm=FALSE)
  if(is.na(f_new[1])) f_new[1] = f_new[2]
  if(figures) {
    plotdata(f_tmp, ft, ftt,
             f_ft, ft_ftt, f_ftt,
             f_ft_bad, ft_ftt_bad, f_ftt_bad, theta)
  }
  mat = matrix(c(f,ft,ftt), ncol=3)
  mat2 = matrix(c(f_ft, f_ftt, ft_ftt), ncol=2, byrow=T)
  list(f = f_new, N = length(makeNA), navalues = makeNA, mat = mat, mat2 = mat2, theta=theta)
}
nk_despike = function(f, just_despike=TRUE, quiet=TRUE) {
  f_mean = mean(f, na.rm=T)
  i=1
  Nold=100
  deltaN = 100
  lambda = sqrt(2*log(length(f)))
  C1 = 1.483
  C2 = 1.483
  f = f - f_mean
  fsigma = MAD(f)
  um = C2 * fsigma * lambda
  good_data = which(f >= -C1*fsigma, f <= C1*fsigma)
  remove = which(abs(f) > um)
  f[remove] = NA
  f[1] = ifelse(is.na(f[1]), 0, f[1]) #初期値がNAの場合エラーを返すのでNAの場合は便宜上0を代入
  f = zoo::na.locf(f, na.rm=FALSE)
  f = f + f_mean
  number_of_spikes = rep(0, 10)
  while(deltaN > 0 & i < 10) {
    f = f - f_mean
    out=despike(f, good_data, figure=F)
    length(out$f)
    dim(out$mat)
    number_of_spikes[i] = out$N
    Nnew = out$N
    f=out$f + f_mean
    # f_mean = mean(f, na.rm=T)
    f_mean = median(f, na.rm=T)
    deltaN = abs(Nold - Nnew)
    Nold = Nnew
    if(!quiet) {cat(paste0(i, ": Number of spikes detected: ", Nnew, "\n"))}
    i = i+1
  }
  spikes_removed = number_of_spikes[1]-number_of_spikes[i-1]
  if(!quiet) {
    cat(paste0("Spikes identified: ", number_of_spikes[1],
               "\nProportion of spikes: ", round((number_of_spikes[1] / length(f)),4)*100,
               "% of data.\n",
               "Spikes removed: ", spikes_removed))}
  if(just_despike) {
    f
  } else {
    list(f=f, out=out)
  }
}

流速データの読み込みとノイズの除去

vector_names = dir("../回流水槽光合成実験5/Data/Velocity/", pattern = "dat", full = T)

cnames_vectrino=c("time", "counter", "status",
                  "xvel", "yvel", "zvel", "z2vel",
                  "amp1", "amp2", "amp3", "amp4",
                  "snr1", "snr2", "snr3", "snr4",
                  "corr1", "corr2", "corr3", "corr4")

vector = lapply(vector_names, read_table , col_names=cnames_vectrino)

vector_test1 = vector("list", length(vector_names))
vector_test2 = vector("list", length(vector_names))
vector_test_result = vector("list", length(vector_names))

for(i in 1:length(vector_names)){
  # Try on a small subset of data
  # Trial run on the first 10000 data points
  vector_test1[[i]] = vector[[i]] %>% slice(1000:11000) #testのため一部抽出
  
  # Remove spikes
  # First clean up the data based on the SNR and CORR
  x = names(vector_test1[[i]] %>% select(dplyr::contains("vel")))
  vector_test1[[i]][,x]=snr_corr_filter(vector_test1[[i]])
  
  # Next run nk_despike() to clean it all up.
  vector_test2[[i]] = vector_test1[[i]] %>% mutate(tau = 1:length(xvel),
                                                   u = zoo::na.locf(xvel,na.rm=F), #NAを直近のNA以外の数値に置き換える
                                                   v = zoo::na.locf(yvel,na.rm=F),
                                                   w = zoo::na.locf(zvel,na.rm=F))
  
  vector_test_result[[i]] = vector_test2[[i]] %>% mutate_at(vars(u,v,w), nk_despike)
}

names(vector_test_result) = paste("", str_extract(vector_names, "[0-9]{1,2}_[0-9]{1,2}_[0-9]{1,3}_[HLM]"), sep = "")

vector_test_result2 = vector_test_result %>%
  bind_rows(.id = "location") %>%
  mutate(location2 = location) %>%
  tidyr::separate(location2, into=c("Hz", "Depth", "Distance", "Density"), sep = '_') %>%
  select(Hz, Depth,
         Distance, Density,
         time,
         xvel, yvel, zvel, z2vel,
         tau, u, v, w) %>%
  mutate(Depth = as.numeric(Depth),
         Hz = as.numeric(Hz)) %>% 
  mutate(Distance = as.numeric(Distance)) %>% 
  arrange(Hz, Distance, Depth, Density) %>%
  tbl_df()

vector_all = vector_test_result2 %>%
  mutate(Hz = ifelse(Hz == 5, 0.5, Hz)) %>% 
  mutate(Density = recode(Density, H = "High", M = "Middle", L = "Low"))

#保存
write.csv(vector_all, "../回流水槽光合成実験5/Modified_Data/vector_all.csv")

流動の統計量の算出

 算出流動に関する統計量は以下とする。

  • 実験条件、測定値ごとの平均値: \(\bar{u}\)
  • 残差: \(u^{\prime}=u-\bar{u}\)
  • 残差の二乗平均: \(\bar{u}^{\prime}=\sum(u^{\prime})^{2}\:/\:n\)
  • TKE(乱流エネルギー): \(TKE=0.5*(\bar{u}^{\prime2}+\bar{v}^{\prime2}+\bar{w}^{\prime2})\)
  • Reynold stress: \(\bar{RS}= -\sum u^{\prime}*w^{\prime}\:/n\)
     
vector_all$Density = factor(vector_all$Density, levels = c("High", "Middle", "Low")) #密度の順位変更

vector_mean = vector_all %>%
  group_by(Hz, Depth, Distance, Density) %>%
  summarise_at(vars(xvel, yvel, zvel, u, v, w),funs(mean, sd)) 

vector_zansa = vector_all %>%
  left_join(vector_mean, by=c("Depth", "Distance", "Hz", "Density")) %>%
  mutate(zansaX = u - u_mean, zansaY = v - v_mean, zansaZ = w - w_mean)

vector_zansa_mean = vector_zansa %>%
  group_by(Depth, Distance, Hz, Density) %>%
  summarise(zansaX_mean = mean(zansaX^2),
            zansaY_mean = mean(zansaY^2),
            zansaZ_mean = mean(zansaZ^2),
            stress = mean(zansaX * zansaZ) * (-1),
            TKE = 0.5*(zansaX_mean + zansaY_mean + zansaZ_mean))

#藻場直前の測定値の鉛直平均を上流流速とする。
upstream_u_mean = vector_mean %>%
  filter(Distance == 26) %>%
  select(Hz, Depth, Density, upstream_u_mean = u_mean) %>% 
  group_by(Hz, Density) %>% 
  summarise(upstream_u_mean = mean(upstream_u_mean))

vector_mean = vector_mean %>% left_join(upstream_u_mean, by = c("Hz", "Density"))
vector_zansa_mean = vector_zansa_mean %>% left_join(upstream_u_mean, by = c("Hz", "Density"))

#測定深度を藻場の中、藻場の境界、藻場の上にグループ分けする。
vector_mean2 = vector_mean %>%
  mutate(Depth2 = ifelse(Depth >= 7 & Depth <= 13, "edge", as.character(Depth)),
         Depth2 = recode(Depth2, "16" = "above", "6" = "in")) 

vector_zansa_mean2 = vector_zansa_mean %>%
  mutate(Depth2 = ifelse(Depth >= 7 & Depth <= 13, "edge", as.character(Depth)),
         Depth2 = recode(Depth2, "16" = "above", "6" = "in")) 

光データと藻場の形状変化データの読み込み及び編集

光データの編集

光(lux)は、ペンダントロガーを2つを、 leading edgeから15 cm, 45 cmの底面に設置して測定した。インターバルは1秒。

#データの読み込み---------------------------------------------------------------------------------------------
fnames = dir("../回流水槽光合成実験5/Data/Light/", pattern = "csv", full = T)
listdata_light = lapply(fnames[-7], read.csv, header = T, skip = 1)
light_ID = read.csv(fnames[7], header = T)
names(listdata_light) = paste("", str_extract(basename(fnames[-7]), "Inoue_[a-z]{2,4}_[0-9]{1,6}_[LMH]"), sep = "")
listdata_light = lapply(listdata_light, function(x){x %>% select(Time = 2, Light_lux = 4)})

#データの編集--------------------------------------------------------------------------------------------------
light = listdata_light %>%
  bind_rows(.id = "filename") %>% 
  tidyr::separate(col = filename, into = c("Inoue", "Position", "a", "Density"), sep = "_") %>%
  select(Position, Time, Light_lux, Density) %>% 
  mutate(Time = as.POSIXct(parse_date_time(Time, orders = "%m/%d/%y %I:%M:%S %p", tz = "JST", locale = "ja_JP.utf8")))

light_ID = light_ID %>% 
  mutate(starttime = as.POSIXct(as.character(Starttime), tz = "JST")) %>% 
  select(-3)

light$Hz = NA
light$Density = NA

removetime = 30 #s 初めの切り取り時間
extracttime = removetime + 20 #s 使用する時間
for (i in 1:nrow(light_ID)) {
  
  light[light$Time >= light_ID$starttime[i] + removetime &
          light$Time <= light_ID$starttime[i] + extracttime, "Hz"] = light_ID$Hz[i]
  
  light[light$Time >= light_ID$starttime[i] + removetime &
          light$Time <= light_ID$starttime[i] + extracttime, "Density"] = paste(light_ID$Density[i])
}

light = light %>%
  na.omit() %>% 
  group_by(Hz, Density, Position) %>% 
  mutate(Elapsed_time = as.integer((Time - min(Time)))) %>% 
  ungroup() %>% 
  mutate(Density = recode(Density, H = "High", M = "Middle", L = "Low"))
         

#平均値、ばらつきを計算-------------------------------------------------------
light_mean = light %>% group_by(Density, Hz) %>% 
  summarise(Light_mean = mean(Light_lux), Light_sd = sd(Light_lux))

藻場の形状変化データの編集

## 海藻の面積データ(cm^2)------------------------------------------
fname = dir("../回流水槽光合成実験5/Data/Area/", pattern = "area_cm.csv", full = T)
Area_cm = lapply(fname, read.csv, header = T)

#データの編集
names(Area_cm) = basename(fname)

Area_cm = Area_cm %>%
  bind_rows(.id = "filenames") %>% 
  tidyr::separate(filenames, into = c("Density", "a", "b"), sep = "_") %>% 
  select(Density, Area_cm = "Area") %>% 
  mutate(Hz = rep(c(0,0.5,1,2,4,6,8,10,15), 3),
         Density = recode(Density, H = "High", M = "Middle", L = "Low")) %>% 
  group_by(Density) %>% 
  mutate(Area_cm_rate = (1-(Area_cm / first(Area_cm)))*100) #藻場面積の減少率、単位は%

データの統合

光合成速度、光、藻場の形状変化、上流流速、藻場の深度グループごとの流速、藻場の高さデータの結合

#藻場(26cm, 41cm)の中(in)、上の淵(edge)、上(above)ごとの流速
Velocity_bed = vector_mean2 %>%
  ungroup() %>% 
  filter(Distance %in% c(41, 71)) %>%
  mutate(Velocity = sqrt((u_mean^2 + v_mean^2 + w_mean^2)/3)) %>%
  group_by(Hz, Density, Depth2, upstream_u_mean) %>%
  summarise(Velocity_bed = mean(Velocity),
            Velocity_bed_sd = sd(Velocity)) 

#藻場(26cm, 41cm)の中(in)、上の淵(edge)、上(above)ごとのTKE
TKE_bed = vector_zansa_mean2 %>%
  ungroup() %>% 
  filter(Distance %in% c(41, 71)) %>%
  group_by(Hz, Density, Depth2, upstream_u_mean) %>%
  summarise(TKE_bed = mean(TKE),
            TKE_bed_sd = sd(TKE)) 

#藻場(26cm, 41cm)の中(in)、上の淵(edge)、上(above)ごとのstress
stress_bed = vector_zansa_mean2 %>%
  ungroup() %>% 
  filter(Distance %in% c(41, 71)) %>%
  group_by(Hz, Density, Depth2, upstream_u_mean) %>%
  summarise(stress_bed = mean(stress),
            stress_bed_sd = sd(stress)) 

all_data = Slope_DO %>% 
  full_join(light_mean, by = c("Hz", "Density")) %>% 
  full_join(Area_cm, by = c("Hz", "Density")) %>% 
  full_join(upstream_u_mean, by = c("Hz", "Density")) %>% 
  full_join(Velocity_bed %>%
              select(-upstream_u_mean) %>% 
              spread(key = Depth2, value = Velocity_bed) %>% 
              rename(Velocity_above = above, Velocity_edge = edge, Velocity_in = 'in'),
            by = c("Hz", "Density")) %>% 
  full_join(stress_bed %>% 
              select(-upstream_u_mean) %>% 
              spread(key = Depth2, value = stress_bed) %>% 
              rename(stress_above = above, stress_edge = edge, stress_in = 'in'),
            by = c("Hz", "Density")) %>% 
  full_join(TKE_bed %>% 
              select(-upstream_u_mean) %>% 
              spread(key = Depth2, value = TKE_bed) %>% 
              rename(TKE_above = above, TKE_edge = edge, TKE_in = 'in'),
            by = c("Hz", "Density"))

#藻場の高さcm
bed_height = read_csv("../回流水槽光合成実験5/Data/Exp_Height_180522.csv") %>% mutate(Height = Height + 1)

all_data = all_data %>% 
  full_join(bed_height, by = c("Hz", "Density")) 

滞留時間と藻場の空隙率(1-solid_volume_fractions)の算出

流速の変化により藻場が変形するため、藻場の体積は変化する。また、密度によって藻場内の海藻の体積も変化する。この2点を滞留時間の計算に用いる。

Length = (90-26)/100 #m 藻場の長さ
Width = 0.30 #m 藻場の幅

all_data = all_data %>% 
  mutate(Total_bed_volume = Height/100 * Length * Width,
         Net_bed_volume = Total_bed_volume - Seaweed_volume_m3,
         Tr_mean = Total_bed_volume / (Velocity_in * Height * Width), #滞留時間の単位はs
         Solid_volume_fractions = Net_bed_volume/Total_bed_volume*100)#空隙率単位は% 

レイノルズ数の算出

レイノルズ数は、流体と物体大きさから、流れの状態を表現する値である。
\(Re=v*L/u\)
\(v: Velocity\)
\(L: 物体の長さ\)
\(u: 動粘度係数\)

#レイノルズ数の計算-------------------------------------
#Re = vL/u v:velocity, L:物体の大きさ
#計算に用いる物体のスケールLは藻場の高さ(m), 海藻間の距離(m)とする。
u = 0.94*10^-6 #m/s 海水25度の時の動粘度係数(Bilger and Atkinson 1992, Kays and Crawford 1993, Thomas and Atkinson 1997)

all_data = all_data %>% 
  mutate(Re_height = ((upstream_u_mean) * (Height))/u,
         Distance_individuals = recode(Density, "High" = 3, "Middle" = 6, "Low" = 8),#個体間の距離cm
         Re_distance = ((upstream_u_mean) * (Distance_individuals/100))/u,
         Re_height_distance = (upstream_u_mean * ((Distance_individuals+Height)/100))/u) %>% 
  return()

#密度の順位変更
all_data$Density = factor(all_data$Density, levels = c("High", "Middle", "Low"))

基本的な結果

藻場の流動環境の結果

藻場のleading edgeを0cmとし、藻場の終わりは、藻場密度がHigh:54cm, Middle:57cm, Low:54cmとした。

#藻場の高さ(流速、密度条件によって変動する)を追加
vector_mean2 = vector_mean2 %>% 
  left_join(bed_height) %>% 
  mutate(`H/h` = Depth/Height,
         Density = factor(Density, levels = c("High", "Middle", "Low")))

vector_zansa_mean2 = vector_zansa_mean2 %>% 
  left_join(bed_height) %>% 
  mutate(`H/h` = Depth/Height,
         Density = factor(Density, levels = c("High", "Middle", "Low")))

#藻場の縁を0cmとする
vector_mean2 = vector_mean2 %>% 
  mutate(Distance_leading_edge = Distance - 28) #leading_edgeは28cm 

vector_zansa_mean2 = vector_zansa_mean2 %>% 
  mutate(Distance_leading_edge = Distance - 28) #leading_edgeは28cm 

#流速条件(Hz)を上流流速の平均値に置き換えた列を作成
Vel_treatment = upstream_u_mean %>% 
  ungroup() %>% 
  group_by(Hz) %>% 
  summarise(Vel_treatment = round(mean(upstream_u_mean), digits = 2))

vector_mean2 = vector_mean2 %>% 
  left_join(Vel_treatment)

vector_zansa_mean2 = vector_zansa_mean2 %>% 
  left_join(Vel_treatment)

流速、密度条件ごとの水路方向流速の時間平均

図上部の数値は、藻場上流流速の平均値、図右のabove, edge, inはそれぞれのキャノピーの上部(16 cm)、縁(変動)、中(6 cm)を示している。
* キャノピーのどの深度においても、藻場上流の流速が増加するにしたがって流速は増加した。
* leading edgeから下流にいくにしたがって、aboveでは流速が増加し、edge, inでは流速は減少する傾向を示した(文献)。
* aboveの上流流速0.39\(m\:/\:sec\)のとき、密度が高いから流速が高くなる傾向示し、下流にいくにつれ、差は大きくなった。
* edge、inの上流流速0.1~0.39\(m\:/\:sec\)のとき、藻場密度が低いとき流速が高くなる傾向をを示した。
* inではleading edgeから下流にいくにしたがって、流速が減少したが、密度が低密度で、減少の仕方が緩やかであった。


流速、密度条件ごとの鉛直方向の流速の時間平均

  • above、edgeでは藻場の直前(x = -2cm)からx = 13 cmで鉛直方向の流速が高くなり、下流にいくに従い、減少する傾向を示した。
  • 上記の傾向は、藻場密度が高いほうが顕著であった。
  • 藻場に海水がぶつかれば鉛直方向の流速が発生するが、aboveの低密度ではあまり顕著ではなく、inでは密度ごとに傾向は異なっていた。


流速、密度条件ごとのレイノルズ応力

レイノルズ応力は、乱流を考慮した流体の運動方程式を立てたときに発生する項であり、渦によって発生する運動量を示している。
* レイノルズ応力は、上流流速が増加するにしたがって、どの密度、深度条件での増加する傾向を示した。
* aboveでは、高密度で高い値を示し、inでは低密度で高い値を示した。


流速、密度条件ごとの乱流エネルギー(TKE)

TKEは、3方向の成分の残渣の二乗平均の和に0.5をかけた数値であり、乱流エネルギーを示す。
* TKEは上流流速が増加するにしたがって、すべての密度条件で増加した。
* above, edgeでは比較的高密度で高い値を示し、inでは、低密度で高い値を示した。


藻場内の平均流速

3方向の流速の二乗平均を示している。
* 密度が高い方が、edge, inにおいて流速が遅くなる傾向がみらた。
* 密度ごとの流速の違いは上流流速が約0.1 m / secから堅調になっていた。


藻場内の平均TKE


藻場内の平均レイノルズ応力


レイノルズ数


滞留時間


その他の要因


藻場の形状変化


空隙率(藻場全体の体積の内、海藻以外の体積の割合)


参考になる図


光合成速度、呼吸速度の結果

呼吸速度と上流流速

 暗条件下における酸素濃度の変動速度から算出した。

呼吸速度は、藻場密度が低いほど大きく、約10 cm/sをピークに低くなる傾向を示した。


純光合成速度と上流流速

 明条件下における酸素濃度の変動速度から算出した。

 純光合成速度は、呼吸速度と同様に約10 cm/sでピークをとり、その後減少した。  しかし、流速が約26 cm/sで全ての密度条件で増加した。純光合成速度は、呼吸速度とは異なり、密度が高い方が、大きい値を示した。


総光合成速度と上流流速

 純光合成速度と呼吸速度の合計値。

GAMによるあてはめ。上流流速と共に総光合成速度が山型に変動している。


各パラメータの相互関係

総光合成速度と他の要因

総光合成速度と滞留時間

藻場内の海水の滞留時間は、藻場全体への向き炭素供給量と負の関係があると考えられ、滞留時間が短くなるほど、光合成速度は増加するはずである。  

GAMによるあてはめ。上流流速と共に総光合成速度が山型に変動している。


総光合成速度と藻場の形状変化

藻場の形状変化は、流速が高くなるにつれて、大きくなった。藻場が折りたたまれることで、光合成速度へどのような影響があるのか?  

GAMによるあてはめ。上流流速と共に総光合成速度が山型に変動している。


Stanによる総光合成速度のあてはめ

本実験では、藻場上流の流速が藻場全体の総光合成速度にどのように影響を与えているかを明らかにするため、応答変数を総光合成速度、説明変数を上流流速として、モデル式へのあてはめを行った。
以下のモデル式を用いた。

 \(Gross\:rate\:=\:b_{0}\:+\:b_{1}*(1-exp(-b_{2}\:/b_{1}\:*\:Speed))\:*\:exp(-b_{3}\:/b_{1}\:*\:Speed)\)
 

Stanのコード


functions{
  vector modelcurve(real b0, real b1, real b2, real b3, vector speed) {
    return (b0 + b1 * (1-exp(-b2 / b1 * speed)) .* exp(-b3 / b1 * speed));
  }
}
data {
  int<lower=0> N;                 // データの数
  int<lower=0> M;                 // 期待値に使う結果の数
  vector[N] speed;                // 観測した速度
  vector[N] gross;                  // 観測した総一次生産量
  vector[4] prior_log_sigma;      // 各パラメータの事前分散
  vector[4] prior_log_mu;         // 各パラメータの事前期待値
}
transformed data {
  vector[N] gross10;
  vector[M] speed_predict;
  real speed_increment;
  gross10 = gross*10;
  speed_increment = max(speed) / M;
  speed_predict[1] = 0;
  
  for(m in 2:M) {
    speed_predict[m] = speed_predict[m-1] + speed_increment;
  }
}
parameters {
  real<lower=0> offset;
  real<lower=0> pmax;
  real<lower=0> alpha;
  real<lower=0> beta;
  
  real<lower=0> sigma;
}
transformed parameters {
  vector[N] fitted;
  fitted = modelcurve(offset, pmax, alpha, beta, speed);
}
model {
  offset ~ lognormal(prior_log_mu[1], prior_log_sigma[1]);
  pmax   ~ lognormal(prior_log_mu[2], prior_log_sigma[2]);
  alpha  ~ lognormal(prior_log_mu[3], prior_log_sigma[3]);
  beta   ~ lognormal(prior_log_mu[4], prior_log_sigma[4]);
  sigma ~ cauchy(0, 2.5);
  gross10 ~ lognormal(log(fitted), sigma);
}
generated quantities {
  vector[M] expected10;
  vector[M] expected;
  vector[N] log_lik;
  expected10 = modelcurve(offset, pmax, alpha, beta, speed_predict);
  expected = expected10/10;
  for(n in 1:N) {
    log_lik[n] = lognormal_lpdf(log(gross10[n]) | log(fitted[n]), sigma);
  }
}

モデル式へのあてはめ

#パッケージの読み込み
library(rstan)

#Stanに渡すデータを準備(Stan側と同じオブジェクト名にする)
Gross_data = all_data %>%
  filter(Hz != 0, type == "GP")

N = nrow(Gross_data)                                      #データの数
speed = Gross_data$upstream_u_mean                        #観測した速度
gross = Gross_data$Slope                                  #観測した総一次生産量 

#密度条件ごとに変化しないパラメーター
M = 50                                                  #期待値に使う結果の数
prior_log_sigma = as.vector(c(1,1,1,1))               #各パラメータの事前分散
prior_log_mu = as.vector(c(1,1,1,1))

#Stanの実行
rng_seed = 1111
nchains = 4
iter = 1000
ncores = 4 

stanout = sampling(nonlinearmodel_gross, 
                   chains = nchains, 
                   iter = iter, 
                   cores = ncores,
                   seed = rng_seed,
                   control = list(adapt_delta = 1))

saveRDS(stanout, file ="../回流水槽光合成実験5/Modified_Data/stanout.obj") #解析結果が保存されたオブジェクトを保存する

Stanの結果

stanout@model_pars
##  [1] "offset"     "pmax"       "alpha"      "beta"       "sigma"     
##  [6] "fitted"     "expected10" "expected"   "log_lik"    "lp__"
pars =  c( "offset", "pmax", "alpha", "beta", "sigma")
pars_fitted =  c( "fitted")
pars_expected =  c( "expected")
traceplot(stanout, pars = c(pars))
## Stan model '8d896e915aa8a88ef1bd2ae55bf0ec2d' does not contain samples.
# traceplot(stanout, pars = c(pars_fitted))
# traceplot(stanout, pars = c(pars_expected))
print(stanout, prob=c(0.025, 0.5, 0.975), par = pars)
## Stan model '8d896e915aa8a88ef1bd2ae55bf0ec2d' does not contain samples.

データの確認